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Correlated electron behavior of metalorganic molecules: insights from density 
functional theory combined with many-body effects using exact diagonalization 
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A proper theoretical description of electronic structure of the 3d orbitals in the metal centers 
of functional metalorganics is a challenging problem. In this letter, we apply density functional 
theory and an exact diagonalization method in a many body approach to study the ground state 
electronic configuration of an iron porphyrin (FeP) molecule. Our study reveals that dynamical 
correlation effects are important, and FeP is a potential candidate for realizing a spin crossover due 
to a subtle balance of crystal field effects, on-site Coulomb repulsion and hybridization between the 
Fe d-orbitals and ligand N p-states. The mechanism of switching between two close lying electronic 
configurations of Fe-d orbitals is shown. We discuss the generality of the suggested approach and 
the possibility to properly describe the electronic structure and related low energy physics of the 
whole class of correlated metal centered organometallic molecules. 


Molecular magnets combine low dimensionality and in¬ 
herent confinement effects with strong electron correla¬ 
tions and hold prospects in the context of spintronics. An 
important molecular property is bistability, i.e., the pos¬ 
sibility of realizing two different spin states, which can in 
principle be accessed and manipulated externally. This 
is important as the switching of spin has a pronounced 
effect on measurable quantities, like magnetic anisotropy 
and spin dipole moment contribution [3[. Finding w^s 
for the efficient manipulation of the magnetic state [l|- 
Q of transition metal (TM) centered porphyrin (TM-P) 
and phthalocyanine (TM-Pc) molecules have critical con¬ 
sequences in this regard. A crucial interplay between 
molecular ligand field and spin pairing energy makes only 
a subspace of this class of materials to respond to spin 
crossover. 

The magnetic properties in TM-P/TM-Pc are largely 
governed by the metal center, which features sizable lo¬ 
cal Coulomb interactions {U ^ 4 eV and J ^ leV) 
and is simultaneously subjected to crystal fields, spin- 
orbit coupling and orbitally dependent hybridization 
with the ligands. Electronic correlations are expected 
to arise 0 and the description by local density approx¬ 
imation (LDA) or generalized gradient approximation 
(GGA) thus potentially becomes inadequate, for exam¬ 
ple, leading to underestimated or even vanishing HOMO- 
LUMO (HOMO=Highest Occupied Molecular Orbital 
and LUMO=Lowest Unoccupied Molecular Orbital) gap. 
Hence, the treatment of the molecular electronic struc¬ 
ture in terms of correlated electron theories like ligand 
field theories or Anderson impurity models become cru¬ 
cial. These model based approaches can be very helpful 
to trace the physical origin of phenomena like spin-state 
switching [HI, emergence of magnetic anisotropies @ or 


many body resonances 0 as soon as solid links between 
model and the realistic structure can be established. 

In this letter, we have adapted a hybrid approach 
0 (DFT++), which links Density Functional Theory 
(DFT) and Anderson’s impurity model to study physical 
properties of FeP and FePc. We demonstrate how the 
interplay of Goulomb interactions, crystal fields and hy¬ 
bridization with the ligands, which are fully captured in 
our theory, lead to correlated electron physics, and how 
this theory describes the S = 1 ^ S = 2 spin-crossover in 
the Fe^^ metal center of FeP. Furthermore, the crossover 
between different close lying ground state electronic con¬ 
figurations within the S = 1 subspace of the correlated 
Fe-d orbitals is analyzed. 

The problem can be cast in the following way. The de¬ 
localized orbitals in the organic ring in FeP are described 
by LDA or GGA but the Fe center is considered as an 
impurity embedded in the organic host and is described 
by Anderson’s impurity model [H. 

The model Hamiltonian for the impurity problem can 
be expressed as : 

H = ^ E Uiiikd\d]dkdi 

i,j hjikd 

+ J2{Vikcldi + H.c..) + J2 ikclck, ( 1 ) 

ik k 

where efj describes the onsite energies and i,j,k,l = 
(m, cr) represents combined orbital and spin indices, di is 
the annihilation operator while Uijik represents the local 
Coulomb interaction. Uijik is parametrized by Slater pa¬ 
rameters. We have chosen U=4 eV and J=1 eV for / = 2 
(3d-orbitals) in our calculations. In the above equation, 
the first two terms represent electrons in Fe 3d orbitals. 
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The third term describes the interaction with the sur¬ 
rounding atoms while the fourth term is for the delocal¬ 
ized ligand electrons with energies The hopping ma¬ 
trix element Vik appears in the hybridization function, 
which is represented as : 


k 


^ik^kj 

e^i5 - ek 


( 2 ) 


The energy dependent hybridization function can be ob¬ 
tained from first principles calculations. We follow the 
approach considered by Karolak et al. 0 to construct 
the local Green function from DFT. The Kohn-Sham 
Green function Gks can be calculated from the Lehmann 
representation [ll| using 


Gks{G) 
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( 3 ) 


where V^n/c’s and e^/c’s are the Kohn-Sham eigenstates 
and eigenvalues for band n and reciprocal space point 
k. The projection of the full Green’s function to an atom 
centred local Green’s function for localized orbitals 

is needed, which in our case are cubic harmonics (Xm)- 


( 4 ) 


where = (xm | ^nk) and P^ = \Xm')- The 

hybridization function is calculated from the local impu¬ 
rity Green’s function from the expression : 

= C + — Vcryst ~ ^(c) = 6 + — A(e). (5) 

In the above expression. Gimp is the projected Green’s 
function on local orbitals and A combines the hybridiza¬ 
tion function A and static crystal field Vcryst • If the bath 
orbitals (defined by ek and Vik in Eqn. ([2])) are limited to 
a small number of discrete orbitals only, the many body 
problem defined in Eqn. m can be solved by means of 
exact diagonalization (ED), which will be used here. 

Electrons in Ee-3d orbitals hybridize mostly with the 
orbitals of the four surrounding N-atoms, that provide 
a square planar ligand field. The effect of the outer G- 
ring is rather indirect as that mainly shifts N-p levels by 
both in-plane and tt — tt interaction. The energy depen¬ 
dence of that interaction with surrounding N-atoms can 
be described by an energy-dependent hybridization func¬ 
tion A(e). We have employed non-spin polarized density 
functional calculations within local density and general¬ 
ized gradient approximations to extract the hybridiza¬ 
tion functions. The DET calculations were performed 
using the VASP code 0 that employs a plane wave ba¬ 
sis and projector augmented wave method. In Eig. [U 
real and imaginary parts of A are shown. The imaginary 


part of ImA quantifies the density of bath states cou¬ 
pling to each impurity orbital weighted by the hybridiza¬ 
tion matrix elements Vik- As seen from Eig. [TJ the most 
dominant peak in ImA(e) is observed for the Ee-da, 2_^2 
orbital at 2.04 eV below the Eermi energy. The forma¬ 
tion of in-plane a bonds with axial N-ligands explains 
this pronounced peak. It should be noted that the other 
in-plane orbital dxy shows almost no hybridization apart 
from a small peak at 4.8 eV below the Eermi energy. The 
out-of-plane orbitals have relatively small peaks in ImA, 
among which the one closest to the Eermi energy is of 
d^^ character at -2 eV. Appearance of this peak reflects a 
TT — TT interaction of Ee d^^ {dxz^ ^yz) orbitals with N-p;^ 
orbitals, which is expected in the square planar ligand 
field of the EeP molecule. The other out-of-plane and 
in-plane contributions are present in the -4.5 to -10 eV 
energy range. Taken together, the hybridization function 
reveals a strong in-plane interaction between Ee-d 2 , 2_^2 
orbitals and the N-p ligand states along with a much 
weaker interaction amongst all other orbitals. 
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FIG. 1. (Color online) Real and imaginary parts of the hy¬ 
bridization function for Fe in FeP calculated with PBE in the 
non-spin polarized mode. The geometry of FeP is shown in 
the inset with the atoms labeled by their types. 


The real part of the hybridization function ReA de¬ 
scribes the energy dependent ligand field, combined with 
a static crystal field, which can be obtained at the e ^ oo 
limit. The strong resonance with the host and da, 2_^2 or¬ 
bital manifests itself also around -2.04 eV in ReA. This 
strong ligand field pushes da, 2_^2 high in energy causing 
almost no occupancy in either spin channel. In the gas 
phase, six electrons in Fe^+, thus, are distributed in the 
remaining four d-orbitals, with four and two electrons in 
majority and minority spin channels, respectively, giving 
rise to an intermediate spin state (S=l). 

An intriguing phenomenon of strain induced spin state 
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switching is observed theoretically in FeP iii , which 
is relatively difficult to obtain for other molecules with 
other TM centers [mil or different structures, e.g., 
transition metal pht halo cyanines. Iron Pht halo cyanine 
(FePc), for example, has an Fe center but along with a 
larger organic ring. This results in a stronger hybridiza¬ 
tion of Fe with the neighboring N atoms, which is re¬ 
flected in the appearance of a dx2_y2 -peak in ImA at 
-1.85 eV and a stronger hybridization, as shown in the 
Supplementary Information (SI). The shift of the bath 
energy to -1.85 eV for FePc compared to -2.04 eV in FeP 
is due to the presence of four extra N atoms connected 
to the pyrolle rings in FePc (See Fig. 2 in SI). These four 
N atoms are not directly bonded to Fe but are connected 
to N atoms in the Fe-N 4 block via C atoms. 

To quantify the conditions required for spin state 
switching of FeP, we have varied crystal field and hy¬ 
bridization strength within DFT+-h method. The varia¬ 
tion of these parameters mimics the strain effect on the 
molecule. Also as mentioned in the previous section we 
found that the crystal field splitting of da, 2_^2 is larger 
compared to other orbitals and as will be discussed be¬ 
low, this splitting is responsible for the spin switching. 
For this part of the calculation we have kept a common 
reference level of the rest of the orbitals separated from 
<1^2 _y2 by Very St ^ which was varied. An independent vari¬ 
ation of the hopping matrix elements, Vd^ 2 _ 2 ? done 
along with the variation of Veryst- One needs to keep in 
mind that for our model calculation, Ycryst and Vd^ 2 _y 2 
are independent, while in DFT calculation (discussed in 
the latter part of this section), these parameters are im¬ 
plicitly related. The ligand field variation due to the 
strain effect simultaneously changes Ycryst and Vd^ 2 _ 2 - 
For the remaining part of the discussion we will refer to 
Vd^ 2 _ 2 as V, as this the only bath-site coupling consid¬ 
ered &r our model calculations. 

Fig. [2] depicts the spin phase diagram of the central 
atom in FeP. The phases are defined by the character¬ 
istic energy contributions and demonstrated with RGB 
color code (see supplementary information). For V 0 
and with sufficiently high crystal field {Ycryst >2.6 eV), 
6 electrons occupy the degenerate drest level. Hence, the 
ground state will be a low spin state with an energy gain 
from crystal field (red) but with the cost of exchange en¬ 
ergy (blue). The ground state attains a high spin state 
for low crystal fields. A strong V pushes the molecu¬ 
lar orbital containing predominantly da, 2_^2 and hence, 
required crystal field for spin crossover is small. The in¬ 
termediate spin state, in this situation is governed by the 
Fe-N bonding (green). Crystal field can be tuned to oc¬ 
cupy anti-bonding molecular orbital and hence, S=2 spin 
state can be achieved (AS). The dependence of V with 
y cry St is follows a quadratic behaviour. The fitted curve 
in Fig. [2] at the phase boundary is obtained from a mean 
field model (see supplementary information) containing 
one particle Hamiltonian but with a renormalized onsite 
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FIG. 2. (Color online) The phase diagram depicting the spin 
states of FeP with the tuning of the static crystal field Ycryst 
and the hybridization strength V. The phase boundary is ac¬ 
companied by the allowed values of Ycryst for fixed values of 
V. The blue curve is a result of fitting with a tight-binding 
model described in the Supplementary Information. AS and 
BS indicate antibonding and bonding regions. Calculated val¬ 
ues of Ycryst and V from DFT (non spin-polarised PBE) are 
shown in orange circles along with the corresponding Fe-N 
bond lengths. 

energy e^. The discrepancy between the phase bound¬ 
ary and the fitted curve lies in proper renormalization of 
e^, that depends linearly with A, by electron correlation 
energy, which has a dependence of N‘^. 

As the spectrum of the FeP molecule is gapped, we 
placed chemical potential in the middle of the gap. But 
it is also possible to vary the average onsite d-energy 
while keeping the total number of particles in the system 
fixed. In addition, the average d-electron on-site energy is 
not exactly known from the DFT simulations — a prob¬ 
lem usually referred to as ’’double counting problem” in 
DFT-h+, DFT+U, or DFT+DMFT approaches. In our 
case, this means that the average on-site energy or more 
precisely the average energy difference between the bath 
level and the Fe d-block carries some uncertainty. We 
thus vary Cd in a range, while keeping the total number 
of particles in our system constant. In this way, the high- 
spin to low-spin transition line has error bars associated 
with some of the calculated points shown in Fig. [2j 

At a stronger coupling V, the shifting of the onsite 
energy ed will allow a variation of Ycryst^ where a spin 
crossover (SGO) can happen. The range of Ycryst due to 
the variation of Cd is presented by the red error bars in 
Fig. [2j In the regime of weak hybridization, the system 
is described predominantly by crystal fields. There are in 
particular, no charge fluctuations to A = 7 or A = 5 Fe 
impurity states, which are generally affected by changing 
ed. As there is no mixing with N ^ d impurity states 
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in the limit of ^ —> 0 , the error bar is vanishingly small 
in this limit. For V beyond 2.8 eV, this particular sce¬ 
nario of crystal field will not be able to switch the spin 
state. An orbital reversal, i.e., the dx 2 _y 2 energy becom¬ 
ing lowered compared to other orbitals is needed in this 
case, which requires a different kind of charge distribu¬ 
tion in the molecule. Thus a transition between S=1 and 
S =2 states could be realized for V < 2.8 eV while S =1 
should be obtained generally for V > 2.8 eV. 
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FIG. 3. (Color online) Magnetic anisotropy energy (MAE) as 
a function of hybridisation V. The corresponding spin states 
of the ground state configurations are indicated, (inset) En¬ 
ergy differences between low spin (LS) and high spin (HS) 
states calculated by DET-h+ and PBE-hU (double counting 
by Dudarev and Lichtenstein) methods as a function of Ee-N 
bond lengths in EeP are shown. The error bars correspond to 
the variation of onsite energy (model Hamiltonian in eq. 1) 
in DET-h+ calculations. 

Weakening of the ligand field leading to a spin state 
change, however, needs a nontrivial chemical or physical 
procedure. The existence of high spin (S= 2 ) porphyrin 
complexes with dP configuration has so far been observed 
in non-planar molecules with five or six coordination of 
the central Fe atom[l5|. Higher coordination leads to out 
of plane shift of the central Fe atom (five coordination) or 
symmetrical Fe-N block expansion, resulting in a weaker 
ligand field. For four coordination, however, the S =2 
state is yet to be observed experimentally in gas phase 
or on a surface. Fig. [2] establishes the parameter space 
for when this state is to be expected. 

The variation of V is also studied with PBE by vary¬ 
ing the Fe-N bond length in the molecule. The high¬ 
est value of V (3.16 eV) is obtained for a Fe-N bond 
length of 1.97 A. The four data points, shown in Fig. [2] 
in filled orange circles, are for bond lengths of 2.00, 2.04, 


2.07 and 2.11 A, respectively. As mentioned in Refs. [H 
and 0 for FeP either physisorbed or chemisorbed on sur¬ 
faces, the required bond length of Fe-N in FeP for spin 
switching is beyond 2.03 A within PBE+U approxima¬ 
tion with Ue//=3 eV. 0 This is in agreement with the 
data shown in Fig.O where both DFT- 1 -+ and PBE sug¬ 
gest a spin crossover beyond 2.04 A. A detailed compari¬ 
son between the bond lengths required for spin-crossover 
in PBE-hU and DFT+-h is shown in Fig. [H A com¬ 
parison of PBE and LDA calculations yield values of 
the static crystal field, bath energy and hybridisation 
as 0.62(0.65) eV, -2.04(-1.9) eV and 3.16(3.39) eV for 
PBE(LDA). These values indicate that the two approxi¬ 
mations for exchange-correlation functional do not show 
pronounced difference. 

In Fig. [3l calculated magnetic anisotropy energies 
(MAE) within DFT-^-h are shown for different strengths 
of V. The MAE is defined as the energy difference be¬ 
tween the lowest many-body eigenstate with out-of-plane 
magnetic moment and the ground state which has in¬ 
plane oriented moment. It is observed that for low val¬ 
ues of V (large Fe-N bond lengths), the MAE decreases 
with decreasing V whereas for large V, a more or less 
constant value of MAE is obtained. In all the cases, 
we have found the easy axis of magnetization to lie in 
the plane of the molecule. In the inset of Fig. 3, the 
spin-crossover properties of FeP as predicted by differ¬ 
ent flavors of PBE-I-U and DFT-l--h is compared. In the 
PBE+U calculations, we stabilize the non-favorable spin 
solutions by constraining the spin moments. It is seen 
from the figure that there is a notable difference between 
two PBE+U methods (Dudarev and Lichtenstein) in the 
Fe-N bond length required for LS-HS transition, while the 
results obtained from Lichtenstein PBE+U and DFT++ 
methods are close to each other. This is understand¬ 
able since both, the Lichtenstein variant of PBE+U and 
DFT+-h, employ the full four fermion Coulomb matrix as 
defined through the Slater parameters, whereas Dudarev 
assumes a simplified Coulomb vertex. 

The free molecule spin state (S=l) is particularly in¬ 
teresting from the point of view of ground state configu¬ 
ration. A strong ligand field in the free molecule leaves 
dx 2 _y 2 nearly unoccupied and an intermediate spin state 
(S=l) can have ground state electronic configuration 
with 2 electrons in dxy, 3 electrons in d^^ and 1 electron in 
dz 2 . We will refer to this configuration a C 231 (d^^d^djs) 
configuration. The other possible configurations within 
the S =1 ground state multiplet are C 222 C 141 

(d^^d^d^s) and C 132 (d^^d^d^s) among which the C 222 
appears to be very close in energy to the energy of C 231 . 
From ReA, a splitting can also be seen among d^^ and 
dxy orbitals but in a relatively small energy scale. To ac¬ 
quire a clear view of the ground state configuration, we 
varied the crystal field in presence and absence of cou¬ 
pling V. At the first step we only considered a crystal field 
splitting between da, 2_^2 and the averaged position of re- 



























5 



FIG. 4. (Color online) Schematic diagram depicting the C222 
and C231 configurations. Ycryst describes the energy separa¬ 
tion between d 3 . 2_^2 and the rest of the d-orbitals treated as 
degenerate. V denotes the strength of hybridization between 
the da, 2_^2 and the bath state with an energy E?,. The specific 
orbitals corresponding to the occupation of the energy levels 
are also indicated. 


maining orbitals, as shown in the left most part of Fig. ID 
In this pure crystal field situation, the S=1 state leads to 
the C 231 ground state configuration. As the coupling V 
is switched on, the ground state configuration becomes 
onsite energy dependent. This happens because the in¬ 
clusion of strong hybridization will admix N=5 and N=7 
impurity occupancies with the pure N =6 ground state 
configuration of the pure crystal field situation. Vary¬ 
ing on site energy, the ground state configuration can 
be modified. The energy scale associated to this change 
is of the order of meV. A relative splitting among the 
remaining orbitals has even more pronounced effects in 
determining the ground state configuration. As shown in 
Fig. m with additional splitting, if d^^ stays above 
C 222 is stabilized. However, in the reversed situation, 
C 231 is obtained. This change occurs in at least one or 
two order higher energy scale compared to the many- 
body effects induced configuration change, revealing the 
crystal field to be the dominant factor. 

In summary, we have presented ground state electronic 
properties of Fe in FeP molecule with a hybrid approach 
of DFT combined with a many body treatment, using 
exact diagonalization. We have demonstrated that a 
delicate interplay of the static crystal field, ligand hy¬ 
bridization and Coulomb interactions promotes iron por¬ 
phyrin to be a potential candidate for realizing spin- 


crossover behavior. In general, our calculated phase di¬ 
agram indicates the possibility of tuning electronic and 
magnetic properties of organometallic molecules to serve 
the purpose of molecular electronics and storage devices. 
Moreover, the long-standing debate regarding the elec¬ 
tronic ground state configuration, e.g., C 222 vs. C 231 
has been solved by identifying the proper parameters re¬ 
quired to switch one to the other, which will have im¬ 
portant consequences for the spin dipole moments and 
magnetic anisotropies, where the energy positions of d- 
orbitals with specific symmetries are important. 
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